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Abstract 

Continuum  damage  formulations  are  commonly  used  for  the  simula¬ 
tion  of  diffuse  fracture  processes.  Implicit  gradient  damage  models  are 
employed  to  avoid  the  spurious  mesh  dependencies  associated  with  lo¬ 
cal  continuum  damage  models.  The  C^-continuity  of  traditional  finite 
elements  has  hindered  the  study  of  higher-order  gradient  damage  approx¬ 
imations.  In  this  contribution  we  use  isogeometric  finite  elements,  which 
allow  for  the  construction  of  higher-order  continuous  basis  functions  on 
complex  domains.  We  study  the  suitability  of  isogeometric  finite  elements 
for  the  discretization  of  higher-order  gradient  damage  approximations. 


1  INTRODUCTION 

Continuum  damage  models  [1]  are  widely  used  for  the  simulation  of  diffuse  frac¬ 
ture  processes.  Several  modifications  of  the  original  theory  have  been  proposed 
to  overcome  the  mesh  dependency  problems  associated  with  the  absence  of  an 
internal  length  scale  (see  e.g.  [2,  3]).  One  way  to  avoid  mesh  dependencies  is 
to  relate  the  material  parameters  to  the  element  size  [4,  5].  Alternatively,  an 
internal  length  scale  can  be  introduced  by  a  spatial  smoothing  function  in  the 
continuum  formulation  [6] .  Gradient  approximations  of  this  smoothing  function 
have  led  to  the  development  of  damage  models  where  an  internal  length  scale  is 
introduced  through  gradients  of  an  equivalent  strain  field.  Among  the  gradient 
damage  formulations,  the  implicit  gradient  enhancement  [7]  is  considered  the 
most  effective.  In  its  original  form  a  second-order  Taylor  expansion  is  used  to 
approximate  a  smoothing  integral,  which  results  in  a  system  of  two  second-order 
partial  differential  equations.  This  formulation  is  attractive  from  a  discretiza¬ 
tion  point  of  view  since  it  can  be  solved  using  C^-continuous  finite  elements. 
It  has,  however,  been  demonstrated  that  the  accuracy  of  the  second-order  ap¬ 
proximation  can  be  limited  [8,  9].  For  that  reason  it  is  important  to  study  the 
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effect  of  the  higher-order  terms  in  the  Taylor  approximation  of  the  nonlocal 
formulation,  which  result  in  higher-order  gradient  damage  formnlations. 

Mixed  finite  element  formulations  can  be  used  for  the  discretization  of  higher- 
order  gradient  damage  formulations.  In  these  formulations,  the  introduction  of 
higher-order  continuous  basis  functions  is  avoided  by  introducing  auxiliary  fields. 
This  results  in  systems  with  many  more  degrees  of  freedom  than  required  by  the 
second-order  gradient  formulation,  making  the  method  computationally  expen¬ 
sive.  To  avoid  the  introduction  of  auxiliary  fields,  meshless  methods  have  been 
used  [9].  The  smoothness  of  meshless  methods  is  inherently  derived  from  the 
way  in  which  the  basis  functions  are  constructed.  Although  meshless  methods 
have  been  applied  successfully  for  the  discretization  of  the  fourth-order  gradi¬ 
ent  damage  formulation,  they  have  not  been  used  widely.  A  reason  for  this  is 
the  inability  of  meshless  methods  to  define  geometry  [10].  The  incompatibility 
with  traditional  finite  element  formulations,  in  the  sense  that  the  method  is  not 
element-based,  may  be  another  reason  why  meshless  methods  are  not  commonly 
applied  to  higher-order  gradient  damage  formulations. 

In  this  contribution  we  use  isogeometric  finite  elements  to  overcome  the 
problems  associated  with  the  use  of  mixed  formulations  and  meshless  meth¬ 
ods  for  gradient  damage  formulations.  The  isogeometric  analysis  concept  was 
introduced  by  Hughes  et  al.  [II]  and  has  been  applied  successfully  to  a  wide 
variety  of  problems  in  solids,  fluids  and  fluid-structure  interactions  (see  [12]  for 
an  overview).  Use  of  higher-order,  smooth  spline  bases  in  isogeometric  analysis 
has  computational  advantages  over  standard  finite  elements,  especially  when 
higher-order  differential  equations  are  considered  [13].  In  contrast  to  meshless 
methods,  the  geometry  and  solution  space  are  one  and  the  same.  This  makes 
it  possible  to  construct  bases  for  complex  geometries,  which  can  be  obtained 
directly  from  a  computer  aided  design  (CAD)  tool  [14].  From  an  analysis  point 
of  view  isogeometric  analysis  can  be  considered  as  an  element-based  discretiza¬ 
tion  technique.  This  compatibility  with  traditional  finite  elements  facilitates 
the  application  to  industrial  problems. 

We  first  review  the  nonlocal  continuum  damage  formulation  and  the  gradient- 
based  approximation  in  section  2.  We  then  introduce  in  section  3  the  isogeomet¬ 
ric  finite  element  discretization.  In  section  4  we  present  numerical  simulations 
utilizing  isogeometric  finite  elements  for  the  discretization  of  the  second-order, 
fourth-order  and  sixth-order  gradient  formulations. 


2  ISOTROPIC  DAMAGE  FORMULATION 


We  consider  a  body  D  C  with  N  G  {1,2,3}  and  boundary  dfl  (see  Figure  1). 
The  displacement  of  a  material  point  x  G  is  denoted  by  u(x)  G  The 
displacements  satisfy  Dirichlet  boundary  conditions,  Ui  =  Ui,  on  dflui  Q  dfl. 
Under  the  assumption  of  small  displacement  gradients,  the  infinitesimal  strain 
tensor 


(1) 


2 


Figure  1:  Solid  domain  fl  with  boundary  dfl. 


is  used  as  an  appropriate  measure  for  the  deformation  of  the  body.  The  Cauchy 
stress  tensor,  a{x)  G  ,  is  used  as  the  corresponding  stress  measure.  An 

external  traction  ti  acts  on  the  Neumann  boundary  dflt^  C  dQ  and  is  equal 
to  the  projection  of  the  stress  tensor  on  the  outward  pointing  normal  vector 
n{x)  €  i.e.  ti  =  o’ij'nj.  The  solid  body  is  loaded  by  increasing  the  boundary 
tractions  or  boundary  displacements. 

2.1  Constitutive  modeling 

In  isotropic  continuum  damage  models,  the  Cauchy  stress  is  related  to  the  in- 
hnitesimal  strain  tensor  by 


—  (1  (2) 

where  w  €  [0, 1]  is  a  scalar  damage  parameter  and  H  is  the  Hookean  elasticity 
tensor  for  undamage  material  (i.e.  with  w  =  0).  When  damage  has  fully  devel¬ 
oped  (w  =  1)  a  material  has  lost  all  stiffness.  Note  that  we  adopt  index  notation 
with  summation  from  1  to  over  repeated  italic  subscript  indices,  for  example. 

The  damage  parameter  is  related  to  a  history  parameter  «;  by  a  monotonically 
increasing  function  lo  =  uj{k),  which  is  referred  to  as  the  damage  law.  Various 
damage  laws  will  be  considered  in  the  numerical  simulations  section.  The  history 
parameter  evolves  according  to  the  Kuhn- Tucker  conditions 

/  <  0,  K  >  0,  k/  =  0  (3) 

for  the  loading  function  f  =  fj  —  k,  where  is  a  nonlocal  strain  measure,  referred 
to  as  the  nonlocal  equivalent  strain.  The  monotonicity  of  both  k  and  uj{k) 
guarantees  that  the  damage  parameter  is  monotonically  increasing  at  every 
material  point,  thereby  introducing  irreversibility  in  the  constitutive  model. 

Nonlocality  is  introduced  into  the  model  by  means  of  the  nonlocal  equivalent 
strain  which  ensures  a  well-posed  formulation  at  the  onset  of  damage  evolution. 
If  instead  the  damage  parameter  was  related  to  a  local  strain  measure,  77,  the 
resulting  medium  would  suffer  from  a  local  loss  of  ellipticity  in  the  case  of 
material  softening  [15].  The  model  is  then  unable  to  smear  out  the  damage 
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zone  over  a  finite  volume.  In  other  words,  a  local  continuum  damage  formulation 
fails  to  introduce  a  length  scale  for  the  damage  zone,  resulting  in  spurious  mesh 
dependencies  in  numerical  solutions. 

A  straightforward  way  of  introducing  nonlocality  in  the  formulation  is  by 
defining  the  nonlocal  equivalent  strain,  fj{x),  as  the  volume  average  of  a  local 
equivalent  strain,  rj  =  77(e), 


r]{x)  = 


f  g{x,y)r]{y)dy 
yen _ 

/  g{x,y)dy 
yen 


where  g{x,y)  is  the  weighting  function 


g[x,  y)  =  exp  - 


\\x-y\Y 

2n 


(4) 


(5) 


We  refer  to  this  model  as  the  nonlocal  damage  formulation  [6] .  The  local  equiv¬ 
alent  strain  maps  the  strain  tensor  onto  a  scalar. 

Although  the  nonlocal  formulation  is  straightforward,  it  requires  the  com¬ 
putation  of  a  volume  integral  for  the  evaluation  of  the  constitutive  behavior 
at  every  material  point.  This  makes  the  numerical  implementation  both  cum¬ 
bersome  and  inefficient.  In  particular,  the  stiffness  matrix  is  full.  Even  when 
truncated,  the  nonlocal  operator  has  a  negative  impact  on  the  sparsity  of  the 
matrix.  This  results  in  computationally  expensive  assembly  and  solution  rou¬ 
tines.  To  circumvent  these  deficiencies,  approximations  of  the  integral  equation 
are  commonly  used. 

The  nonlocal  equivalent  strain  (4)  can  be  approximated  by  substitution  of  a 
Taylor  expansion  for  the  equivalent  strain  field  around  the  point  x 


{y^-x^){yj-XJ)  +  0{{xi-yif). 

y^x 

(6) 

Assuming  the  solid  volume  stretches  to  infinity  leads  to  the  gradient  approxi¬ 
mation  of  equation  (4) 


v{y)  =  ^ly=x+^ 


ri{x)  =  r]{x)  +  -I, 


2 "  ax? ' 


dS 

dx'ldx'] 


+  48^' 


a®  77 


’ dxldx^dxl 


(x)  +  ....  (7) 


This  gradient  approximation  is  known  as  the  explicit  gradient  formulation.  As 
an  alternative,  the  implicit  gradient  formulation  (e.g.  Ref.  [7])  is  obtained  by 
direct  manipulation  of  equation  (7) 


n(ct 


1  .n  a^n 


Because  only  C^-continuity  is  required  for  the  second-order  approximation,  the 
corresponding  implicit  gradient  formulation  has  enjoyed  widespread  use. 
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In  the  remainder  of  this  work  we  study  the  convergence  of  the  implicit  gra¬ 
dient  formulation  toward  the  nonlocal  formulation  upon  increasing  the  number 
of  gradient  terms  involved.  If  we  truncate  equation  (8)  after  the  d-th  derivative, 
we  can  rewrite  it  using  a  linear  operator  as 

C‘^fi{x)  =  ri{x).  (9) 

We  restrict  ourselves  to  the  second-order  {d  =  2),  fourth-order  (d  =  4)  and 
sixth-order  (d  =  6)  implicit  gradient  damage  formulations. 

2.2  Implicit  gradient  damage  formulation 

In  contrast  to  the  nonlocal  and  explicit  gradient  damage  formulations,  the  im¬ 
plicit  formulation  requires  the  solution  of  a  boundary  value  problem  for  the 
nonlocal  equivalent  strain  field,  fj^x),  in  addition  to  the  usual  problem  for  the 
displacement  held,  u{x).  In  the  absence  of  body  forces,  the  resulting  boundary 
value  problem  for  the  d-th  order  formulation  is  given  by 


d(Jii 

n 

dxj 

Vx 

G 

=  V 

dijUj 

=  U 

Vx 

G 

d^u 

-9-  ( 

dXn  \ 

d'^fj  ^ 
dxj...  J 

> 

0 

II 

G 

90,  a  G  {0, . . . ,  d  —  2} 

= 

Ui 

Vx 

G 

90„, 

where  t  and  u  are  the  prescribed  boundary  traction  and  displacements,  respec¬ 
tively.  Notice  that  we  assume  all  directional  derivatives,  of  the 

nonlocal  equivalent  strain  held  are  zero  on  the  boundary.  We  verify  this  choice 
numerically  by  comparing  the  results  with  the  nonlocal  formulation  based  on 
the  integral  equation  (4).  The  kinematic  and  constitutive  relations  (1)  and  (2) 
are  used  to  express  the  Cauchy  stress  in  terms  of  the  displacement  held. 

We  solve  the  system  (10)  using  the  Galerkin  method.  The  same  solution 
spaces  are  used  for  the  displacement  held  and  nonlocal  equivalent  strain  held, 
denoted  by  5“  C  and  C  id  2(0)^  respectively.  We  denote  our  trial 

spaces  as  V“  and  and  assume  that  and  V“  and  5“  are  the  same 

modulo  inhomogeneous  boundary  conditions.  The  weak  form  of  equation  (10) 
then  follows  as 


V?;“  e  V“ 

Vu’'  G 


(11) 


where  ^  and  (•,  •)n  is  the  L^-inner  product.  No  boundary 

terms  appear  in  the  equation  for  the  equivalent  strain  held,  since  the  deriva¬ 
tives  of  this  held  in  the  direction  of  the  normal  vector  are  assumed  zero  on  the 
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boundary  of  the  domain.  For  the  damage  formulations  considered  in  this  work 
(i.e.,  with  d  G  {2, 4,  6}),  the  linear  operator  is  written  as 


JilA 

^2  dxi’ 


^dxidxf  y/^dxidxjdxk 


(12) 


3  ISOGEOMETRIC  FINITE  ELEMENTS 

Discretization  of  the  weak  formulation  (11)  for  the  d-th  order  damage  formu¬ 
lation  requires  (|  —  l)-times  continuously  differentiable  basis  functions.  With 
isogeometric  finite  elements  -continuous  basis  functions  can  be  constructed 
using  splines  of  order  p.  We  will  consider  B-splines  and  non-uniform  rational 
B-splines  (NURBS)  [16].  Although  T-splines  [17]  are  considered  beyond  the 
scope  of  this  contribution,  it  is  emphasized  that  a  unihed  analysis  approach  for 
splines  is  provided  by  Bezier  extraction  [18]. 


3.1  Univariate  B-splines  and  NURBS 

The  fundamental  building  block  of  isogeometric  analysis  is  the  univariate  B- 
spline,  e.g.  [16,  12].  A  univariate  B-spline  is  a  piecewise  polynomial  defined  over 
a  knot  vector  S  =  {^i,  ^2,  •  •  ■ ,  Cn-i-p-i-i},  with  n  and  p  denoting  the  number  and 
order  of  basis  functions,  respectively.  The  knot  values  are  non-decreasing 
with  increasing  knot  index  i,  i.e.  Ci  ^  '?2  <  •  ■  •  <  Cn-i-p-i-i-  A  partition  of  the 
parameter  space  [^i,^„+p+i)  is  provided  by  the  elements 

Se  =  [?*(e),'Ci(e)-Hl)  (13) 


where  z(e)  is  a  map  from  the  element  indices  to  the  knot  vector  indices. 

The  B-spline  basis  is  defined  recursively,  starting  with  the  zeroth  order  [p  = 
0)  functions 


1 

0  otherwise 


(14) 


from  which  the  higher-order  [p  =  1,2,...)  basis  functions  can  be  constructed 
using  the  Cox-de  Boor  recursion  formula  [19,  20] 


W,p-i(g)  +  ^  jV.+i.p-i(e)-  (15) 

S*-l-p  si  S*-l-p-l-l  S»-l-l 

Efficient  and  robust  algorithms  exist  for  the  evaluation  of  these  non-negative 
basis  functions  and  their  derivatives,  e.g.  [21].  B-spline  basis  functions  satisfy 
the  partition  of  unity  property  and  B-spline  parameterizations  possess  the  vari¬ 
ation  diminishing  property,  e.g.  [22].  An  example  of  a  univariate  B-spline  basis 
is  shown  in  Figure  2.  For  notational  convenience,  we  will  drop  the  subscript  p 
of  the  basis  functions. 
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A  B-spline  is  defined  as  a  linear  combination  of  B-spline  basis  functions 

n 

=  (16) 

i=l 

where  Ai  is  called  a  control  point  or  variable.  Equation  (16)  is  typically  used  for 
the  parameterization  of  curves  in  two  (with  €  R^)  or  three  (with  Xi  G  K^) 
dimensions  by 

n 

*(e)  =  (17) 

i=l 

In  combination  with  equation  (13),  this  parametric  map  provides  a  definition  of 
the  elements  in  the  physical  space 

o,  =  {x(e),  C  e  sj .  (18) 


B-splines  used  for  analysis  purposes  are  generally  open  B-splines,  i.e.,  the  first 
and  last  knot  values  have  a  multiplicity  of  p+  1.  Dirichlet  boundary  conditions 
can  be  applied  to  the  control  points  on  the  boundary  of  the  domain  in  the  same 
way  as  done  for  the  nodal  variables  in  traditional  finite  elements.  B-splines 
can  also  be  refined,  which  is  important  in  the  context  of  isogeometric  analysis, 
e.g.  [23]. 

A  drawback  of  B-splines  is  their  inability  to  exactly  represent  many  objects 
of  engineering  interest,  such  as  conic  sections.  For  this  reason  NURBS,  which 
are  a  rational  generalization  of  B-splines,  are  commonly  used.  A  NURBS  is 
defined  as 

n 

«(e)=E7?*(^)^-  (19) 

with  the  NURBS  basis  functions  defined  as 


MO 


MQWp 


(20) 


where  w{^)  =  Ni{^)Wi  is  called  the  weighting  function.  Note  that  in 

equation  (20)  no  summation  is  performed  over  the  repeated  index  j3.  In  the 
special  case  that  Wi  =  c  Vi  G  {l,...,n},  and  c  an  arbitrary  constant,  the 
NURBS  basis  reduces  to  the  B-spline  basis. 


3.2  Multivariate  B-splines  and  NURBS 

Computational  domains  in  two  and  three  dimensions  can  be  parametrized  by 
means  of  bivariate  and  trivariate  splines,  respectively.  An  A^-variate  B-spline  is 
contracted  over  the  tensor  product  of  knot  vectors 

(g)  . . .  0  =  1^]^, . . . , }  0  ...  0  , . . .  .  (21) 
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Figure  2:  Third  order  B-spline  basis  for  the  global  knot  vector  S  = 
{0,0, 0,0, 1,2,  3, 4, 4, 4, 4}. 


where  Ua  and  Pa  are  the  number  of  basis  functions  and  order  of  the  basis 
functions  in  the  direction  a,  respectively.  The  parameter  space  0 

■  ■  ■  ®  [^i  ,^nN+PN+i)  partitioned  by  the  elements 


where  ia{e)  is  a  map  from  the  element  indices  to  the  indices  of  the  knot  vector 
2“  for  a  =  1 ... 

Multivariate  B-spline  basis  functions  with  parametric  coordinate  ^  = 

. . .  ,C^),  are  defined  through  the  product 


N 

m)  =  n 

a=l 


(23) 


where  jaii)  maps  the  global  control  point  indices  on  the  indices  of  the  univariate 
basis  functions  in  the  direction  a.  Multivariate  NURBS  basis  functions  are 
defined  by  equation  (20)  using  a  weights  vector  W  G  M",  where  n  =  0^=1 
is  the  number  of  multivariate  basis  functions  . 

Using  the  multivariate  basis  functions,  an  A"-variate  B-spline  is  defined  by 
equation  (19).  When,  in  addition  to  a  weights  vector  W,  a  set  of  control  points 
X  G  is  provided,  a  parameterization  of  an  Wdimensional  body  11  is 

obtained  by 

n 

3^(0  =  (24) 

i=l 

In  combination  with  the  partitioning  of  the  multivariate  parameter  space  by 
the  elements  2e  in  equation  (22)  this  parametric  map  provides  the  definition  of 
elements  He  in  tbe  physical  space  through  equation  (18). 


3.3  Isogeometric  finite  element  discretization 

Let  5“^  c  5“  and  C  be  the  discrete  solution  spaces  for  the  displace¬ 
ment  held,  u{x),  and  nonlocal  equivalent  strain  held,  fj{x),  respectively.  In 
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isogeometric  analysis,  these  spaces  are  written  in  terms  of  the  NURBS  basis 
functions  Ri{^)  defined  in  section  3.1  and  3.2  for  the  univariate  and  multivari¬ 
ate  cases,  respectively.  We  can  approximate  the  displacement  held  and  nonlocal 
equivalent  strain  held  as 


=  '^Rk{^{x))Uki 

k^l 

n 

=  '^Rk{i{x))Hk 

fc=l 


(25) 


where  U  G  are  the  control  point  displacements,  and  H  G  M"  are  the  con¬ 

trol  point  nonlocal  equivalent  strains.  From  the  displacement  held,  the  strain, 
e{x),  and  local  equivalent  strain,  r](x),  can  be  computed.  In  combination  with 
the  nonlocal  equivalent  strain  held,  the  damage  parameter,  uj(x),  and  Cauchy 
stress,  cr(x),  can  be  obtained  at  every  point  using  the  constitutive  relations 
provided  in  section  2.1. 

We  use  the  Galerkin  method  to  discretize  the  weak  formulation  (11)  as 


G  v; 


(u.h  \  (7  u.h  \ 

=  Joo 


(26) 


Using  the  NURBS  basis  functions,  Ri{^),  as  trial  functions  results  in  a  system 
of  (TV  -I-  l)n  equations 


/i“nu  =  /exqfe  V(A:,m)e{l...n}®{l...TV} 
fL,k  =  0  VfcG{l...n} 


which  can  be  solved  for  every  load  step  using  Newton-Raphson  iteration  to  de¬ 
termine  the  control  point  coefficients  Uki  and  Hk  in  equation  (25).  The  internal 
force  vectors  are  assembled  by  looping  over  the  elements 


Tie 

rUm  _  A 


= 

J  int,A; 


A 


int.fe 


A 

e=l 


1  fdRk 

Cij ,  _ 


2  \  dxj 

rie 


^in 

d/2 

E 

13=1 


dRk. 

dx, 


(28) 


The  integrals  in  these  expressions  are  evaluated  on  the  elements  Ug  dehned 
through  equation  (18).  In  this  contribution,  we  use  Gaussian  quadrature  of 
order  p  -|-  1  in  each  direction.  Numerical  integration  of  NURBS  for  analysis 
purposes  was  studied  in  [24]  and  remains  an  active  topic  of  research.  In  order 
to  evaluate  the  integrals,  the  Jacobian  of  the  isogeometric  map  needs  to  be 
evaluated  at  every  integration  point.  Since  rational  basis  functions  are  used. 
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Figure  3:  Schematic  representation  of  a  one-dimensional  rod  loaded  in  tension. 
The  cross-sectional  area  of  the  rod  is  10  mm^  except  for  the  central  section  where 
it  is  equal  to  9  mm^ . 


this  requires  application  of  the  quotient  rule.  Since  higher-order  derivatives 
with  respect  to  the  physical  coordinate  x  are  used  in  this  contribution,  higher- 
order  derivatives  of  the  parametric  map  are  also  required.  These  higher-order 
derivatives  are  obtained  by  application  of  the  chain  rule  to  the  differentiation 
of  the  NURBS  basis  functions  with  respect  to  the  physical  coordinates  [12]. 

The  consistent  tangent  matrix,  required  by  the  Newton-Raphson  procedure, 
can  be  obtained  by  differentiation  of  (28)  with  respect  to  the  control  point 
variables  in  equation  (25).  Evaluation  of  the  tangent  requires  the  derivatives 
of  the  stress,  tr,  with  respect  to  the  strain,  e,  and  nonlocal  equivalent  strain,  fj. 
These  derivatives  are  provided  through  the  constitutive  behavior  elaborated  in 
section  2.1.  The  derivative  of  the  local  equivalent  strain,  rj,  with  respect  to  the 
strain  tensor  follows  from  the  equivalent  strain  law,  77  =  r](e). 


4  NUMERICAL  SIMULATIONS 


4.1  One-dimensional  rod  loaded  in  tension 


We  consider  a  one-dimensional  rod  loaded  in  tension  as  shown  in  Figure  3.  The 
10  mm  wide  central  section  of  the  rod  has  a  reduced  cross-sectional  area  in  order 
to  develop  a  centralized  damage  zone.  The  modulus  of  elasticity  of  the  rod  is 
E  =  20  GPa,  and  the  Cauchy  stress  is  written  as  cr  =  (1  —  uj)Ee.  As  a  damage 
law  we  consider  [7] 


uj{k) 


I  0 


K  <  Kq 


K  —  Kq 
Ku  —  I^Q 


K>  Ko 


(29) 


with  Ko  =  1  •  10“^  and  k„  =  0.0125.  We  define  the  local  equivalent  strain  law 
as  ry  =  (e)  where  (•)  is  the  Macaulay  bracket  and  take  the  nonlocal  length  scale 
in  (5)  equal  to  C  =  V^mm. 

Force-displacement  curves  have  been  determined  for  the  second-,  fourth-  and 
sixth-order  implicit  gradient  models  and  for  the  nonlocal  damage  formulation.  A 
dissipation-based  path-following  constraint  [25]  is  used  to  trace  the  equilibrium 
path  beyond  the  snapback  point.  Mesh  convergence  studies  have  been  performed 
using  uniform  meshes  with  80,  160,  320,  640  and  1280  elements  of  orders  1,  2 
and  3.  The  control  points  are  equidistantly  spaced  and  all  control  weights  are 
equal  to  1,  which  leads  to  a  linear  parameterization  of  the  domain.  An  overview 
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of  the  meshes  is  given  in  Table  1.  Note  that  in  contrast  to  higher-order  finite 
elements,  the  number  of  degrees  of  freedom  is  practically  independent  of  the 
order  of  the  basis  (since  for  practical  meshes  p  -C  rie). 

The  force-displacement  results  for  the  second-order  gradient  formulation  ob¬ 
tained  on  all  meshes  are  presented  in  Figure  4.  Increasing  the  order  of  the  basis 
from  linear  to  quadratic  significantly  improves  the  rate  of  convergence.  A  further 
increase  of  the  order  of  the  basis  to  cubics  has  minor  effects  on  the  mesh  conver¬ 
gence  behavior.  This  is  explained  by  the  fact  that  both  quadratic  and  cubic  ba¬ 
sis  functions  are  unable  to  accurately  represent  the  solution  in  the  region  where 
loss  of  ellipticity  has  occured.  Accurate  results  for  the  second-order  formulation 
are  obtained  using  1280  elements  of  any  order.  Force-displacement  results  for 
the  fourth-order  gradient  formulation  are  shown  in  Figure  5.  The  presence  of 
fourth-order  spatial  derivatives  in  this  formulation  requires  -continuity.  For 
that  reason,  meaningful  results  are  only  obtained  on  quadratic  and  cubic  meshes. 
The  response  obtained  on  the  1280  element  meshes  cannot  be  visually  distin¬ 
guished  from  that  obtained  on  the  640  element  meshes.  A  very  accurate  result 
is  already  obtained  on  the  320  element  mesh.  The  improved  convergence  behav¬ 
ior  of  the  fourth-order  formulation  compared  to  the  second-order  formulation  is 
attributed  to  the  fact  that  smoother  results  are  obtained.  The  increased  smooth¬ 
ness  of  the  fourth-order  formulations  compared  to  the  second-order  formulation 
is  closely  related  to  the  postponed  loss  of  ellipticity  for  these  formulations  as 
demonstrated  by  the  dispersion  analysis  performed  in  [9] .  Meaningful  results  for 
the  sixth-order  formulation  are  only  obtained  on  cubic  meshes,  and  are  shown 
in  Figure  6.  The  convergence  behavior  of  the  sixth-order  formulation  closely 
resembles  that  of  the  fourth-order  formulation.  In  Figure  7  we  present  the  re¬ 
sults  obtained  on  all  meshes  for  the  nonlocal  formulation.  Compared  to  the 
results  obtained  using  the  second-order  formulation,  the  nonlocal  formulation 
shows  superior  convergence  behavior  for  the  linear,  quadratic  and  cubic  meshes. 
As  expected,  the  difference  in  convergence  behavior  diminishes  as  the  order  of 
the  gradient  formulation  is  increased.  In  fact,  the  results  for  the  nonlocal  for¬ 
mulation  obtained  using  cubic  basis  functions  can  hardly  be  distinguised  from 
those  obtained  using  the  sixth-order  formulation  (Figure  6).  For  the  purpose 
of  comparing  the  various  formulations,  the  accuracy  of  the  solutions  obtained 
on  the  1280  cubic  element  meshes  suffices.  A  detailed  study  of  the  convergence 
rates  is  a  topic  of  future  research. 

In  Figure  8  we  show  a  comparison  of  the  various  formulations.  All  results 
are  obtained  on  a  mesh  with  1280  cubic  elements.  The  results  are  in  excellent 
agreement  with  those  reported  in  e.g.  [7]  and  [9].  As  in  [9]  it  is  observed  that 
the  incorporation  of  fourth-order  derivatives  in  the  implicit  scheme  improves 
the  results,  in  the  sense  that  the  obtained  force-displacement  curve  is  closer  to 
that  of  the  nonlocal  formulation.  Consistent  with  this  observation  we  find  that 
the  sixth-order  formulation  gives  an  even  better  approximation  of  the  nonlocal 
result.  Similar  trends  are  observed  from  the  final  damage  profiles  (see  Figure  8). 
The  sixth-order  formulation  is  found  to  be  very  efficient  since  the  results  are 
in  good  agreement  with  the  nonlocal  formulation,  while  the  involved  computa¬ 
tional  effort  is  very  small  compared  to  the  nonlocal  formulation.  Based  on  the 
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Linear  B-spline  (p  =  1) 


Number  of  elements,  Ug  80 

160 

320 

640 

1280 

Number  of  basis  functions,  n  81 

161 

321 

641 

1281 

Quadratic  B-spline  (p  = 

2) 

Number  of  elements,  Ug  80 

160 

320 

640 

1280 

Number  of  basis  functions,  n  82 

162 

322 

642 

1282 

Cubic  B-spline  (p 

=  3) 

Number  of  elements,  Ug  80 

160 

320 

640 

1280 

Number  of  basis  functions,  n  83 

163 

323 

643 

1283 

Table  1:  Meshes  used  for  the  uniaxial  rod  simulation.  Note  that  for  this  case 
n  =  He  +  P- 


(a)  (b) 


(c) 


Figure  4:  Mesh  convergence  studies  for  the  rod  using  the  second-order  gradient 
formulation,  discretized  using  linear  B-splines  (a),  quadratic  B-splines  (b),  and 
cubic  B-splines  (c).  The  key  labels  indicate  the  number  of  elements. 
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u  [mm] 


u  [mm] 


(a) 


(b) 


Figure  5:  Mesh  convergence  studies  for  the  rod  using  the  fourth-order  gradient 
formulation,  discretized  using  quadratic  B-splines  (a),  and  cubic  B-splines  (b). 
The  key  labels  indicate  the  number  of  elements. 


0  0.02  0.04  0.06  0.08 

u  [mm] 


Figure  6:  Mesh  convergence  studies  for  the  rod  using  the  sixth-order  gradi¬ 
ent  formulation,  discretized  using  cubic  B-splines.  The  key  labels  indicate  the 
number  of  elements. 
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(a)  (b) 


(c) 


Figure  7:  Mesh  convergence  studies  for  the  rod  using  the  nonlocal  formulation, 
discretized  using  linear  B-splines  (a),  quadratic  B-splines  (b),  and  cubic  B- 
splines  (c).  The  key  labels  indicate  the  number  of  elements. 
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u  [mm] 


Figure  8:  Force-displacement  diagrams  for  the  rod  loaded  in  tension  using  the 
nonlocal  formulation  and  d-th  order  gradient  formulations.  All  results  are  ob¬ 
tained  using  1280  cubic  elements. 


Figure  9:  Final  damage  profile  for  the  rod  loaded  in  tension  using  the  nonlocal 
formulation  and  the  d-th  order  gradient  formulations.  All  results  are  obtained 
using  1280  cubic  elements. 


resemblance  of  the  sixth-order  and  nonlocal  result  it  is  concluded  that,  for  the 
considered  simulation,  ignoring  the  nonlocal  equivalent  strain  boundary  terms 
appearing  in  the  gradient  formulation  has  a  minor  effect  on  the  results. 

4.2  Three-point  bending  beam 

As  a  second  numerical  experiment  we  consider  the  three-point  bending  beam 
experiment  considered  in  [9]  and  shown  in  Figure  10.  The  1000  x  300  mm^  beam 
is  supported  by  hinges  on  the  left  and  right  bottom  corners,  and  is  loaded  by  a 
distributed  load  t  over  the  central  100  mm  section  of  the  specimen. 

A  linear  isotropic  material  is  considered  with  modulus  of  elasticity  E  = 
20GPa  in  the  undamaged  state  and  Poisson’s  ratio  ly  =  0.2.  Plane  strain  con- 
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Figure  10:  Three-point  bending  specimen.  The  thickness  of  the  specimen  is 
50  mm. 


ditions  are  assumed,  and  the  local  equivalent  strain  is  given  by 

vis)  =  {Sif  =  \/  {sif  +  {£2f  (30) 

where  ei  and  £2  are  the  principal  strains  of  the  two-dimensional  strain  tensor. 
The  Macaulay  bracket  distinguishes  the  cases  of  tension  and  compression.  The 
following  damage  law  as  proposed  in  [26]  is  used 


0  K  <  kq 

1  {(1  “  “)  + [/^(^o  -  k)]}  k>ko 


(31) 


with  parameters  kq  =  1  •  10“^,  a  =  0.99  and  /3  =  500.  The  nonlocal  length  scale 
is  taken  as  Ic  =  20  mm. 

Force-displacement  curves  are  obtained  using  the  meshes  shown  in  Figure  11 
with  linear,  quadratic  and  cubic  basis  functions.  A  summary  of  the  mesh  pa¬ 
rameters  is  given  in  Table  2.  Note  that  in  two  dimensions  the  number  of  basis 
functions  is  also  practically  independent  of  the  order  of  the  formulation,  in  con¬ 
trast  to  traditional  cubic  finite  elements.  The  control  point  weights  are  all  taken 
equal  to  1,  and  control  points  are  placed  such  that  a  bilinear  parameterization  of 
the  beam  is  obtained.  The  force  F  is  defined  as  the  distributed  load  t  times  the 
width  to  which  it  is  applied.  The  displacement  u  is  measured  by  the  downward 
displacement  of  point  A,  as  indicated  in  Figure  10.  This  displacement  is  used 
as  a  path-following  constraint  to  trace  the  equilibrium  path. 

In  Figure  12  we  show  the  response  curves  for  the  second-order  gradient  for¬ 
mulation  obtained  on  all  meshes.  As  in  the  one-dimensional  study  in  section  4.1 
we  see  a  signihcant  improvement  in  the  convergence  rate  when  using  quadratic 
basis  functions  instead  of  linear  basis  functions.  Accurate  results  are  obtained 
using  Mesh  4  with  quadratic  or  cubic  basis  functions.  The  results  obtained  on 
the  quadratic  and  cubic  meshes  for  the  fourth-order  formulation  are  shown  in 
Figure  ??.  For  both  the  quadratic  and  cubic  basis  functions,  the  results  ob¬ 
tained  on  Mesh  4  cannot  be  visually  distinguished  from  those  obtained  on  Mesh 
3.  The  same  holds  for  the  results  obtained  for  the  sixth-order  formulation,  as 
shown  in  Figure  14.  Due  to  the  considerable  computational  effort  required,  the 
mesh  convergence  behavior  of  the  nonlocal  formulation  is  only  studied  using 
Mesh  1,  2  and  3.  As  seen  from  Figure  15,  a  significant  change  in  the  force- 
displacement  curves  is  observed  when  comparing  the  Mesh  2  and  Mesh  3  results 
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Mesh  1 


Mesh  3 


Figure  11:  Meshes  for  the  three-point  bending  specimen. 


Linear  B-spline  (p  =  1) 


Mesh  1 

Mesh  2 

Mesh  3 

Mesh  4 

Number  of  elements,  rig 

128 

512 

2048 

8192 

Number  of  basis  functions,  n 

153 

561 

2145 

8385 

Quadratic  B-spline  (p  =  2) 

Mesh  1 

Mesh  2 

Mesh  3 

Mesh  4 

Number  of  elements,  rig 

128 

512 

2048 

8192 

Number  of  basis  functions,  n 

180 

612 

2244 

8580 

Cubic  B-spline  (p 

=  3) 

Mesh  1 

Mesh  2 

Mesh  3 

Mesh  4 

Number  of  elements,  rie 

128 

512 

2048 

8192 

Number  of  basis  functions,  n 

209 

665 

2345 

8777 

Table  2:  Meshes  used  for  the  three-point  bending  beam.  Note  that  for  these 
meshes  n  =  {\/2n^  +  -|- p). 
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(c) 


Figure  12:  Mesh  convergence  studies  for  the  three-point  bending  beam  using 
the  second-order  gradient  formulation,  discretized  using  linear  B-splines  (a), 
quadratic  B-splines  (b),  and  cubic  B-splines  (c)  . 


with  linear  basis  functions.  When  using  quadratic  basis  function,  the  response 
curves  cannot  be  visually  distinguised.  For  the  purpose  of  comparing  the  various 
formulations,  the  result  of  the  nonlocal  formulation  obtained  on  the  quadratic 
Mesh  3  is  sufficiently  accurate. 

In  Figure  16  the  results  of  the  various  formulations  are  compared.  Typical 
damage  profiles  are  shown  in  Figure  17.  Upon  increasing  the  order  of  the  for¬ 
mulation  the  approximation  of  the  nonlocal  result  is  improved.  Increasing  the 
order  of  the  formulation  increases  the  total  amount  of  dissipated  energy.  This 
is  caused  by  the  additional  smoothing  effect  of  the  higher-order  derivatives.  For 
the  considered  problem,  the  sixth-order  formulation  is  observed  to  be  very  effi¬ 
cient,  since  it  accurately  approximates  the  nonlocal  result,  whereas  the  involved 
computational  effort  is  negligible  compared  to  the  nonlocal  formulation.  As  in 
the  case  of  the  rod  simulation,  setting  all  the  Neumann  boundary  conditions 
(10)  for  the  equivalent  strain  field  to  zero  does  not  have  a  significant  effect  on 
the  results. 
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(a) 


(b) 


Figure  13:  Mesh  convergence  studies  for  the  three-point  bending  beam  using 
the  fourth-order  gradient  formulation,  discretized  using  quadratic  B-splines  (a), 
and  cubic  B-splines  (b). 

labelfig:threepointconvergence4th 


Figure  14:  Mesh  convergence  studies  for  the  three-point  bending  beam  using 
the  sixth-order  gradient  formulation,  discretized  using  cubic  B-splines. 


(a) 


(b) 


Figure  15:  Mesh  convergence  studies  for  the  three-point  bending  beam  using 
the  nonlocal  formulation,  discretized  using  linear  B-splines  (a),  and  quadratic 
B-splines  (b). 
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Figure  16:  Force-displacement  results  for  the  three-point  bending  beam  speci¬ 
men  using  the  nonlocal  formulation  and  d-th  order  gradient  formulations.  All 
results  for  the  gradient  formulation  are  obtained  using  Mesh  4  with  cubic  basis 
functions.  The  result  for  the  nonlocal  formulation  is  obtained  on  Mesh  3  with 
quadratic  basis  functions. 


Figure  17:  Damage  profiles  for  the  three-point  bending  experiment  obtained 
using  the  fourth-order  gradient  formulation.  Undamaged  material  is  indicated 
in  blue,  fully  damaged  material  in  red.  Deformations  are  50  times  amplified. 
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5  CONCLUSIONS 


Isogeometric  analysis  allows  for  the  construction  of  smooth  basis  functions  on 
complex  domains,  providing  an  appropriate  solution  space  for  higher-order  dif¬ 
ferential  equations.  Dirichlet  boundary  conditions  can  be  applied  by  specifying 
control  variables  along  the  boundary,  in  the  same  way  as  nodal  variables  are 
specified  for  traditional  finite  elements. 

Isogeometric  analysis  is  shown  to  be  a  very  good  candidate  for  the  dis¬ 
cretization  of  higher-order  gradient  damage  formulations.  Using  cubic  basis 
functions  allows  for  the  discretization  of  the  sixth-order  gradient  damage  for¬ 
mulation.  Since,  from  a  practical  point  of  view,  the  number  of  degrees  of  freedom 
is  independent  of  the  polynomial  order  of  the  basis  functions,  the  fourth-  and 
sixth-order  formulation  require  only  slightly  more  computational  effort  than  the 
second-order  formulation.  This  makes  it  practical  to  study  the  convergence  of 
the  implicit  gradient  damage  formulation  toward  the  nonlocal  formulation  upon 
increasing  its  order. 

Numerical  simulations  have  been  performed  for  a  one-dimensional  rod  loaded 
in  tension,  for  which  a  univariate  B-spline  basis  is  constructed.  A  two-dimensional 
three-point  bending  beam  specimen  is  discretized  using  a  bivariate  NURBS 
patch.  For  both  simulations  it  is  observed  that  the  result  of  the  nonlocal  for¬ 
mulation  is  approached  upon  increasing  the  order  of  the  gradient  damage  for¬ 
mulation.  Since  the  computational  effort  involved  in  the  nonlocal  formulation 
is  much  larger  than  that  for  the  gradient  approximations,  increasing  the  order 
of  the  gradient  formulation  yields  efhcient  approximations  of  the  nonlocal  re¬ 
sult.  For  the  two  simulations  considered,  the  sixth-order  formulation  gave  an 
accurate  approximation. 
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